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Abstract 

The dynamics of complex systems often involve thermally activated barrier crossing events that allow 
these systems to move from one basin of attraction on the high dimensional energy surface to another. Such 
events are ubiquitous, but challenging to simulate using conventional simulation tools, such as molecular 
dynamics. Recently, Weinan E et al. [Nonlinearity, 24(6), 1831(201 1)] proposed a set of dynamic equations, 
the gentlest ascent dynamics (GAD), to describe the escape of a system from a basin of attraction and proved 
that solutions of GAD converge to index- 1 saddle points of the underlying energy. In this paper, we extend 
GAD to enable finite temperature simulations in which the system hops between different saddle points on 
the energy surface. An effective strategy to use GAD to sample an ensemble of low barrier saddle points 
located in the vicinity of a locally stable configuration on the high dimensional energy surface is proposed. 
The utility of the method is demonstrated by studying the low barrier saddle points associated with point 
defect activity on a surface. This is done for two representative systems, namely, (a) a surface vacancy and 
ad- atom pair and (b) a heptamer island on the (111) surface of copper. 



1 



I. MOTIVATION 



The dynamics of many complex systems proceed via sequences of infrequent transition events 
from one metastable state to another. Well-known examples include chemical reactions, confor- 
mation changes of bio-molecules, nucleation events during phase transition, etc.^ ^ Metastability 
is characterized by the appearance of disparate time scales. Two most important time scales are 
the relaxation time scale of a state, and the transition time out of the state. A state is metastable 
if the second time scale is much bigger than the first. For this reason, transitions between differ- 
ent metastable states are rare events. If the system has a rather simple potential energy landscape 
(PES), the bottleneck for the transition is a saddle point (generally of index-1), called the transition 
state, of the potential energy of the system. 

Crippen and Scheraga pioneered numerical algorithms for climbing out of basins of attraction.l^l 
In their method, starting from an initial point xq, in the k-ih step, the system translates along the 
direction r = (x/^ — xq) by a suitable step length e to yield x^ = x/^ + er. This is followed by 
minimization on a hyperplane Sk perpendicular to r to yield x^^+i = argmin V (x^). This process 
is repeated until the system reaches a saddle point. To sample different saddle points, the authors 
suggested using the eigenvectors {wi, W2, W3, wat} of the Hessian at xq, i.e. Xq = xq + ew^, 
i G [1, A^], followed by minimization on So to obtain xi. 

Later works on exploring the high dimensional potential energy surface (PES) are based on 
the realization that evaluating minimum eigenmode of the Hessian is central to the convergence 
to an index-1 saddle point. The system can however, start from the initial locally stable point 
by following any eigenmode of the Hessian or other randomly selected direction vectors.'^ For 
smaller systems, an eigenvalue problem of the Hessian is easy to solve by diagonalization (though 
one has to remove the global rotational and translational components), but this becomes computa- 
tionally challenging for systems of higher dimensions. Consequently, many techniques have been 
proposed to evaluate the minimum eigenmode of the Hessian in an optimal fashion. 

For instance, Cerjan and Miller suggested a technique that requires selecting a trust region 
around a point on the multidimensional PES and approximating the energy of the system within 
this trust region by a quadratic expression.l^ An optimal direction to translate the system is then 
determined by evaluating the extremum of the energy on the boundary of the trust region. The 
key to proper evaluation of the optimal direction lies on proper selection of the trust region - if 
the trust region too small, then the linear term dominates and only one minimum is obtained; on 
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the other hand, for a large trust region the quadratic approximation of the PES becomes poor. The 
computational effort, for this method, scales up rapidly as more degrees of freedom are included 
in the system and at the same time the construction of the Hessian can become a formidable 
computation challenge when analytical second derivatives are not available. 

A more economic approach to evaluate the ascent direction, that avoids the constructing the 
full Hessian, is the dimer method proposed by Henkelman and Jonsson.*^ A dimer consists of 
two images of the configuration separated by a small distance 10"^ A). The dimer translates 
towards a saddle point by the modified force Fr — 2FII, where Fr is the net force acting on the 
dimer and F" is the component of the force parallel to a dimer direction n. The direction n is 
determined by minimizing the dimer energy, i.e. the direction vector is updated after k-ih iteration 
by 

rifc+i = argmin [V (x/^ + en) + V (x/e - en)] (1) 

n, \n\^=l 

where, V (x) is defines the PES. In the above, V (x/e + en) and V {-kj. — en) represents the energy 
of the dimer separated by a distance of 2e. The process of obtaining the minimum eigenmode 
of the Hessian is implemented by first selecting a suitable plane spanned by the rotational force 
acting on the dimer and the direction n. The dimer is then rotated in this plane to obtain an optimal 
n corresponding to a minimum energy of the dimer.l^ A dynamical version of the dimer method 
was recently proposed by Poddey et al.^ 

Similarly, Munro and Wales suggested an iterative scheme to calculate the minimum eigen 
mode of the Hessian.*^ Starting from an initial guess direction no, the eigenvector corresponding 
to the minimum eigenvalue of the Hessian is found approximately by successive operation of the 
shifted Hessian i.e. n/^ = (H — AI)^ no, where, A is a constant and I is the identity matrix. In 
this procedure, since we only need the product of the Hessian with a vector, there is no need to 
explicitly calculate the Hessian thus reducing the computational complexity significantly. Once 
the minimum eigenmode is determined the system is then translated in the configuration space by 
a suitable step length. 

Other PES exploration algorithms like the activation-relaxation method (ART) proposed by 
Barkema and MousseaiP do not guarantee visiting an index- 1 saddle point. Analysis of the con- 
vergence of ART to saddle points is detailed in the Appendix. Even though ART is an excellent 
exploration tool, the saddle points sampled may not correspond to the low barrier saddle points that 
drive rare transition events. Proposals for fixing this problem can be found in recent literature.'^^HIll 
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The philosophy that underiies the gentlest ascent dynamics (GAD) is to take one step back 
and formulate a set of dynamic equations whose solutions converge to saddle points. Numerical 
algorithms can then be constructed by discretizing this set of dynamic equations. As we illustrate 
in this paper, having a set of dynamical equations helps us to analyze the stability and convergence 
close to different fixed points on the PES. At the same time, it is easy to extend such a formalism 
to finite temperature cases. In this paper, we will discuss how GAD can be naturally extended as a 
tool for sampling saddle points and for the exploration of the configuration space. In particular, we 
will propose finite temperature molecular dynamics version of GAD, MD-GAD. MD-GAD can be 
used to produce trajectories that hop from saddle point regions to saddle point regions. Such time 
series can be useful for a variety of purposes. We also discuss practical issues related to GAD such 
as selecting initial conditions, etc. for high dimensional systems. 

The arrangement of the paper is as follows : in Section II, we discuss the set of dynamical 
equations corresponding to finite temperature and zero temperature conditions to sample saddle 
points. In Section III, we discuss the different ways to implement the set of dynamical equations, 
namely, identifying an important region associated with a rare transition event, proper sampling 
of direction to move out of the energy well and selection of convergence criteria. Next, we use 
two examples -{i) a heptamer island on (111) surface of copper and (n) a surface vacancy and 
ad-atom pair on the (111) surface of copper to study the relevant saddle configurations and the 
energy barriers. 

II. EQUATIONS OF MOTION 

We assume that the system being considered is defined by an energy function V (x) : 
R, where is the number of atoms in the system with a total of 3A^ degrees of freedom (DOF). The 
potential function is assumed to be smooth. The force at a point x G R^^ is F (x) = — (x) 
and the Hessian is given by H (x) = (x). The equations that govern the gentlest ascent 
dynamics are as follows: 

X = F (x) - 2 (F (x) , n) n 

(2) 

7n = — Hn + (n, Hn) n 
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The first equation means that we reverse the force in the direction of n: F = — Fy, where 
F^ and Fy, respectively, are the components of the local force perpendicular to n and parallel to 
n and are given by Fy = (F, n) n, F^ = F — Fy. The second equation defines the dynamics 
of the ascent direction n. The first term on the right hand side ensures that n converges to an 
eigenvector associated with the smallest eigenvalue of H. The second term ensures that the length 
of n is fixed at 1. Note that n (t = 0) should be a unit vector. If H is fixed, then the second 
equation can be regarded as a continuous version of the power method for finding the eigenvector 
corresponding to the smallest eigenvalue. 7 is a parameter that determines the rate of convergence 
of the direction vector to the lowest eigen-mode of the Hessian. In the limit 7^0, the system 
evolves by following the minimum eigenmode of the local Hessian. 

It was shown by E et al that the only stable fixed points of GAD are the index- 1 saddle points.f^ 
In many ways, GAD is the analog of the steepest decent dynamics: 

x = F(x) (3) 

While the steepest decent dynamics are stuck at the local minima, GAD gets trapped at the 
index- 1 saddle points. For this reason, it is of interest to have molecular dynamics or Langevin 
dynamics versions of GAD, which will allow us sample to saddle point regions. This is the one of 
main purposes of the current paper. 

The equations of motion for the stochastic version of GAD are simply: 

X = F (x) - 2 (F (x) , n) n + e r/ (t) 

(4) 

7n = — Hn + (n, Hn) n 

where, r] (t) is the noise term with variance given by (ry (t) ry (t^)) = 5 {t — f) and mean value 
{v (t)) = 0. Since a system following ^ gets trapped at an index- 1 saddle points, presence of a 
noise term in ([4]) helps the system to hop between different saddle points. 

The dynamical system in ([2]) can be extended to molecular dynamics at finite temperature con- 
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ditions as: 



X = V 

V = m-^ [F (x) - 2 (F (x) , n) n] (5) 
7n = — Hn + (n, Hn) n 

Similar to GAD, for the set of dynamical equations ([5]), index- 1 saddle points of V (x) become 
linearly stable and locally stable fixed points of V (x) become linearly unstable fixed points of 
([5]). A detailed mathematical analysis of the local convergence is presented in Appendix I. The 
dynamical system defined by ([5]) is similar to molecular dynamics, where the system hops between 
different basins of attractions around metastable fixed points of V (x), except that in ^ the basins 
of attractions are index- 1 saddle points of V (x). The molecular dynamics version of GAD an also 
be extended to perform constant temperature dynamics: 

X = V 

V = m-^ [F (x) - 2 (F (x) , n) n] - (6) 
7n = — Hn + (n, Hn) n 

where, /5 is a temperature controller similar to Nose-Hover thermostat used in molecular 
dynamicsP'^ 

It is often desirable to perform simulations under constant external stress conditions and obtain 
the ensuing free energy barrier corresponding to a rare event. To this end, we present a variation 
of MD-GAD based on Parrinello-Rahman extension of molecular dynamics to simulate a system 
under constant external stress S and external pressure Let hg be the matrix formed by the 

three initial lattice vectors a, b and c, i.e. hg = {a, b, c} and f^o = det (hg). Correspondingly, 
the position vectors x can be represented in terms of reduced coordinates x^ by the transformation 
x^ = h"^x, where h is the time dependent metric tensor. The equations of motion are 

= h-^ [F - 2 (F, n) n] - Q-^Gv^ 

(7) 

W^h = (H - a - hS 
7n = — Hn + (n, Hn) n 
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In the above, a = dVt/d\\, G = h'^h and W determines the rate at which the system approaches 
equihbrium between the external and internal stress!^ Descriptions for S and 11 can be found in 
Ref[19]. 

A. Illustrative 2-dimensional example 

We start with a 2-dimensional toy example with the energy function : 

V (x, y) = sin (ttx) sin {iiy) , x,y e [-1, 1]. (8) 

V (x, y) has many locally stable points, such as (— |, |), (|, — |), etc. and saddle points (0, 0), 
(±1, 0), etc. inside the domain of interest. First set of simulations are performed using the de- 
terministic version of GAD prescribed in ([2]). The simulation is started at a point close to the 
metastable fixed point (|, — |) with a randomly initialized direction vector. Fig. [l] shows the evo- 
lution of the system in the two dimensional configuration space. As the system moves out of the 
energy well, the direction vector slowly relaxes to the smallest eigenvalue of the Hessian and then 
converges to the saddle point at (-1, 0). An important conclusion that can be drawn from these 
results is that the path can deviate significantly from the minimum energy path and the choice of 
the initial direction vector determines the location of the converged saddle point. Consequently 
different realizations starting from the same initial point can converge to different saddle points. 




FIG. 1 : Potential energy surface for the 2-dimensional toy problem showing the evolution (solid blue line) 
of the system from a metastable state to an index- 1 saddle point by following ([2]). 

The stochastic version of GAD can be used to navigate the PES to sample different saddle 
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points. Fig. |2 (a) [ shows the results for stochastic dynamics performed by the system on the PES. 
The system spends majority of its time around saddle points in contrast to general overdamped 
Langevin dynamics in which the system samples stable fixed points of the PES. 

Next, we turn to MD-GAD prescribed in ([5]). In Fig. 2(6) the system travels from the initial 
locally stable fixed point at (|, — |) to the final metastable state at (|, — |) by passing through 
many saddle points and locally stable fixed points. The simulation is started with zero kinetic 
energy. As the system walks out of the initial energy well and moves towards the saddle point, 
the KE of the system increases and attains a maximum value. This pushes the system towards the 
next metastable point where the KE reduces to zero. Subsequently, the system again gains kinetic 
energy as it moves close to the next saddle point. As this process is repeated many fixed points of 
the PES are sampled. 




o ® o © o © o 
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(a) 



(b) 



FIG. 2: |2(a) Solutions of ^ showing the trajectory followed by a system on the 2-dimensional PES. The 
system spends majority of its time near saddle points. |2(6)| Shows the evolution of the system using MD- 
GAD prescribed in ([5]). In this case the system can easily traverse through many locally stable and index- 1 
saddle points of V (x, y). 



III. APPLICATIONS 

To further illustrate our approach, we simulate an ensemble of saddle points for two systems: 
(a) a heptamer island on (111) surface of copper and (b) a surface vacancy and ad-atom pair on the 
(111) surface of copper. The heptamer island system has been widely studied by Henkelman et 
al. using the dimer method and most recently by Jonsson et al. to benchmark the convergence of 
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different saddle point finding strategies The second example has clear time-scale separation 
between different rare events, namely, ad-atom diffusion and vacancy diffusion and hence presents 
excellent opportunity to test the scope of saddle point finding algorithms. While diffusion of an 
ad-atom on a metal surface has been studied before the presence of a vacancy adds significant 
complexity to the problemF^ Our intention is to explore the high dimensional PES close to these 
locally stable configurations to probe and catalogue some of the probable surface activity that can 
play an important role in other scenarios where quantitative analysis of different rare events may 
be difficult. 

A. Selecting an "active region" of the sample 

The problem of determining saddle points close to a local minimum entails moving out of the 
potential energy well along preferential directions in the multidimensional configuration space. 
The challenge in sampling these preferential directions comes from the fact that they represent 
only a miniscule fraction of all the possible directions emanating from the initial locally stable 
fixed point on the high dimensional PES. 

Our aim is to significantly reduce the size of the configuration space by selecting a small, yet 
important, region of the sample. In materials physics, many deformation processes can be easily 
traced by analyzing the softest eigenmodes of the Hessian. This helps in determining a vulnerable 
region of the system. In systems where one has an intuitive understanding of the vulnerable region 
(Qy) one can select a set of atoms and include their first nearest neighbors to form a set of "active 
region" such that C ^. We show that a more prudent way to sample Qy, which is prone to 
deformation, is by analyzing the unit force vector close to the initial locally stable configuration. 
Once the vulnerable DOE are identified, the corresponding atoms and their first nearest neighbors 
are selected to form the "active region" of the sample. Our assumption is that this "active region" 
corresponds to a significantly reduced configuration space where probable low barrier transition 
events can take place. 

Let us consider the system containing a surface vacancy and an ad-atom on the (111) surface 
of a copper thin film (Eig. |3]). Eor this system, we have explicitly evaluated Hessian for the relaxed 
initial configuration (force norm < 10~^ eV/A). The elements in the i-th row of the Hessian matrix 
are determined by displacing the i-ih DOE by a small displacement Sxi from the equilibrium 
structure and evaluating the change in the force on all atoms. The Hessian (characterized by the 
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upper inset in Fig. 4(a)| ) is then diagonalized and its eigenvalues calculated. Fig. 4(a) shows 
the spectrum of the Hessian of the system. The lower inset in the figure shows the clustering 
of low-lying eigenvalues which may be difficult to resolve by using power method based tools. 
A superposition of the low-lying eigenvectors (first 6 lowest eigenvectors, neglecting the global 
rotation and global translation components) are plotted in Fig. 4(6)[ 

Fig. 4(6)1 also shows the contributions of the different DOF to the unit force vector (shown 
in blue) at the initial locally stable fixed point of the energy landscape. Upon comparison with 
the superposition of the six low-lying eigenvectors, we find the the DOFs that have significantly 
higher weights in the unit force vector also correspond to those in the low lying eigenvectors of the 
Hessian. This suggest that close to a locally stable fixed point the force can provide information 
about the vulnerable DOF of the system and can be used to obtain probable initial direction vectors 
for probing saddle points on the PES. Following the above argument, using the unit normed force 
at the local minimum point, we selected a set of 23 atoms as our "active region" Q. 

We find that even for complicated non-equilibrium processes involving many potential transi- 
tion events in presence of external fields (such as stress, etc.) the unit force on the system can 



be very informative. To illustrate this we use the example of nanoindentation. Fig. 5(a) shows 
a representative system containing 112,320 copper atoms. The simulation cell has lattice vectors 
along [110], [111] and [112] directions. While free surface conditions are prevalent along [111] di- 
rection, periodic boundary conditions are imposed along the perpendicular directions. The system 
is elastically indented using a smooth spherical indenter (radius 25 A), imposed by a force field, to 
an indentation depth of 4.19 A. After indentation, constrained minimization (force norm < 10~^ 



eV/A) is performed to obtain a locally stable configuration. Fig. 5(6) shows the contributions of 
the different DOF in the unit force vector from which we can easily select Analysis of the 
probable transition events taking place during nanoindentation will be presented elsewhere. 



B. Initial direction vectors 



Since the overall system size corresponds to a large configuration space, an activated region 
of the sample is selected following a procedure described in the previous section. From the set 
of atoms within we obtain a set of initial direction vectors by the following procedure : select 
m atoms and randomly initializing the 3m direction components while keeping the remaining 
components of the direction vector equal to zero. The set of m atoms is selected by picking an 
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atom i and its (m — 1) nearest neighbors from Q. If the atom i has n (< m — 1) nearest neighbors 
inside Q then the direction has only (3n + 3) randomly initialized components and the remaining 
DOF are set to zero. For our study of rare events in the representative systems described earlier 
we use m = 1, 2, 4, 12, though one can resort to a different set of values. 

Our choice of selectively activating only a few components of the direction vector stems from 
the fact that even though physical processes, such as different deformation mechanisms, taking 
place in most materials have activation volumes in the range of 0.16^-10006^, 6 is the Burgers 
vector of the system, the size of the nucleus in the initial stages can be quite small. A smaller 
activation volume corresponds to a localized process involving only few atoms such as interstitial 
diffusion, vacancy diffusion etc., while a larger activation volume corresponds to a process which 
is more delocalized, such as, forest hardening, Orowan loopin, ttc^ In each of these processes, 
since only a few degrees of freedom of the direction vector are explored, ideally one needs humon- 
gous number of randomly initialized vectors to sample majority of the most probable rare events. 
Our proposed selection procedure overcomes this sampling problem. 

C. Convergence criteria 

Once a set of the desired direction vectors is obtained, one can move out of basins of attraction 
by moving along these directions. GAD guarantees that the system converges to an index- 1 saddle 
point on the energy surface. However, the convergence will depend on how fast the system can 
relax to the lowest eigenmode in the vicinity of the saddle point. During this time it is possible that 
the system might have traversed through few different fixed points on the PES. In order to avoid 
this "overshooting" one might resort to a smaller time step. However, this is not conducive as the 
majority of initialized direction vectors have sufficient time to relax to the minimum eigenmode 
much before the system reaches the separatrix. 

Hence, a pragmatic approach is needed to determine whether the system has left the basin of 
attraction or not. Our approach is motivated by the existence of an unique equilibrium bond length 
(To) between nearest neighbor atoms in a locally stable configuration. A rare event involves a 
process of bond-breaking between nearest neighbors (which corresponds to a bond length, r) in 
the group of atoms taking part in the process. This criteria can be used to identify the approach of 
a system to a saddle point. We use the ratio q = {r — To) /to for this purpose. To sample saddle 
points close to a local minimum, suggested values of qc lie in the range ^ 0.30 — 0.40.'^^ Using 
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this criteria has an added flexibihty : one can selectively tune qc to determine a distribution of 
saddle points within a given distance from the initial configuration in the hyper-space. 

For the sampling saddle points in the representative examples we use qc = 0.50. The system 
follows ([2]) with 7 = 1 for g < To identify any bond breaking event, q is calculated every 10 
steps. Once, q > qcWt use 7 = 0.25 and for each evolution step of x, the direction vector n is 
evolved 2 times to converge to the minimum eigenmode. 



D. Surface vacancy ad-atom pair on (111) surface of cooper 



The simulation setup consists of a, 480 Cu atoms, simulation cell with axis vectors along [110], 
[111] and [112] directions. The set up consists of 6 (111) planes stacked on top of each other 
with periodic boundary conditions imposed along [110] and [112] directions while free surface 
conditions maintained along [ill] direction. On atom is removed from the surface to form a 
surface vacancy and placed on a vacant fee site close to the vacancy to form an ad-atom. Fig. 
[3] shows the initial configuration obtained by local minimization to a force of 10"^ eV/A. The 
ad-atom is shown in grey. Fig. 4(6)| shows the weights of different DOF in the unit normed force 
vector for the initial configuration. The vulnerable DOF are easily identifiable. These DOF are 
used to generate 661 initial direction vectors. 

Fig. 6(a) to Fig. |7((j)| shows the probable transition events with increasing activation energy 
barriers. The smallest barrier (0.01 eV) corresponds to an ad-atom diffusion. The difference in 
energy between a fee and hep site on a free surface is much smaller than that in metals like Pt, 
resulting in a significantly smaller ad-atom diffusion barrier. We found the process of diffusion of 
both ad-atom and vacancy in a concerted way, as shown in Fig. 6(6)| has a little higher barrier of 
0.15 eV. The other significant events are ad-atom diffusion to annihilate the ad- atom and surface 



vacancy pair (Fig. 6(c)), vacancy diffusion on the surface (Fig. 6{d)) and sub-surface vacancy 
formation ( Fig. 6(e) and |6(/)| ). Fig. 7(a) corresponds to ad-atom migration by exchange process. 
Some of the transition events are affected by the location of ad-atom and vacancy. For example, 
as shown in Fig. 7(6) [ and Fig. |7(c)[ the formation of a divacancy and a nearest- neighbor ad-atom 
pair has lower barrier than formation of a divacancy and two ad-atoms. Fig. |7((j)| corresponds to 
formation of a bulk vacancy and an ad-atom. 

Fig. [8] shows a summary of the results obtained from converged GAD simulations. The results 
are depicted as the distribution of the converged saddle energies as a function of distance of the 
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converged saddle configuration from the initial configuration in the configuration space. As can 
be seen from the figure, there are multiple saddle configurations with same energy but separated 
by a finite distance in the configuration space. 

Fig. |9] shows the distribution of converged saddle energies less than 2.0 eV for different values 
of qc for a fixed set of initial direction vectors. This suggests that, for a sufficiently large set of well 
sampled direction vectors, the values of qc in the range of 0.30-0.50 provides a good collection of 
low barrier saddle points. Though our selection of initial direction vectors significantly decreases 
repeated visits of the same saddle point, however on a high dimensional PES there exists saddle 
points with similar barriers but separated by a finite distance. The above results have significant 
contributions from such scenarios. 



Fig. 10 shows a collection of saddle points obtained by using MD-GAD. In this case, the 



simulations are started with zero initial velocity from the local minimum in Fig. |10(a) and the 



direction vector preferentially initialized to activate the ad-atom. As the simulation progresses, the 
system moves from the local minimum point to a saddle point (Fig. |10(b) ) where the KE reaches 



a maximum value. Then the high KE pushes the system to the next saddle point shown in Fig. 



10(c) Thus MD-GAD effectively samples many low lying saddle points on the PES. 



E. Heptamer island on (111) surface of copper 

The simulation setup consists of a sample with 487 Cu atoms and lattice vectors parallel to 
[110], [ill] and [Il2] directions. The periodic boundary conditions are imposed along [110] and 
[Il2] while free surface conditions are maintained along [ill] direction. The atomic interactions 
are modeled using an EAM potential developed by Mishin et al.f^ Initially all atoms in the seven 
atom heptamer island occupy fee sites on the (111) surface. 

Fig. ll(a)| to Fig. 12((j)| shows a collection of, low barrier, saddle configurations. From our 



simulations, we find many collective processes, such as those involving lateral translation (energy 
barrier 0.39 eV) and rotation (energy barrier 0.92 eV) of the heptamer island. These collective pro- 
cesses have also been reported for islands of different sizes and shapes.^^ Other low barrier events 
include sliding of two nearest neighbor atoms belonging to the heptamer island (Figs. |11(6)[ ll((j)[ 



11(e) I and 12(c)[ ) all with barriers 0.49-0.69 eV. Fig. |ll(/)| involves a cooperative rearrangement 



and sliding of three atoms and has little higher barrier. The sliding of group of atoms with respect 
to others in the island resembles the process of formation of a dislocation in a bulk crystal. We 
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also observe movement of an island by repeated shearing of different layers!^ Fig. Il3| and Il2(rf) 



shows different ways of formation of a surface vacancy and an ad-atom. While the process in 
Fig. 1 2 ((j) I involves a surface atom hopping out of plane to attach itself to the island the process in 
Figjis] involves a cooperative "exchange" mechanism involving a sub-surface atom.f^ 



Fig. [14] shows a summary of the converged saddle configuration energies, obtained using the 
843 initial direction vectors, as a function of the distance in the configuration space from the 
initial local minimum configuration. Setting the parameter = 0.50 provides the flexibility 
to sample a larger portion of the configuration space near the initial configuration. This also 
allows us to capture events which involve multiple saddle points thus providing vital information 



about possible channels for probable transition events. Fig. [15] shows one such example where a 
sequence of transition events leads to the formation of a bulk vacancy, an ad-atom on top of the 
heptamer island and subsequent change in shape of the island on the (111) surface. 



IV. CONCLUSIONS 



We have shown that the deterministic, molecular dynamics or stochastic versions of GAD can 
be used as effective tools for sampling saddle points on the PES. The algorithms presented here 
are scalable to systems of higher dimensions. The computational cost and memory requirements 
are similar to molecular dynamics simulations. 

A method to initialize the direction vectors along which the system can move out of the initial 
potential energy well is described. It is based on the knowledge of the vulnerable DOF of the sys- 
tem from the force at initial local minimum. All the operations in GAD and its finite temperature 
variants have O {N) (N is the number of DOF) computational complexity. Our use of a criteria 
based on equilibrium bond length of atoms in a system to identify the closeness of a system from 
the initial local minimum point also saves computational cost. 

The finite temperature versions of GAD open avenues for sampling of saddle points on the PES. 
Generally existing macroscopic models are based on information obtained from locally stable fixed 
points of a system. In the stochastic version of GAD, since the system can hop between different 
low energy saddle points, a wealth of information can be obtain about the saddle configurations. 
This information can then be used to improve the existing mean field models. 

GAD can be extended to handle higher index saddle points, as was done in Ref[15]. It is 
obvious that one can also extend the work presented to those cases. 
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Our results for point defect activity on a (111) surface obtained from GAD show a variety of 
processes involving more than a single atom. Indeed the wide spectrum of diffusion processes 
revealed by GAD for ad-atoms, vacancies and islands can provide much needed insight into the 
probable mechanisms behind many experimentally observed surface phenomenal^SEll 



V. APPENDIX I:LOCAL CONVERGENCE ANALYSIS FOR MD-GAD 

The fixed points of PES are also the fixed points of ([5]). To see this consider a two dimensional 
system with force F = f^i + fyj and direction vector n = Uxi + Uyj, where fx, fy and Ux, 
riy are individual components of force and direction vectors, respectively, along the i and j axis. 
Consequently, the solutions of v = are: 



fxi + fyj = 2A ( rixi + nyj) , A = /^n^ + fyUy , + = 1 



(9) 



Equating the individual vector components on both sides of the above equation we obtain 



nf, — n: 



:) =2/, 

;) 



yTiyTix 



f! + fy=0 



(10) 



which is possible only if individual components of the force vector are zero. So for all values of 
the direction vector the above condition is satisfied only at fixed points of V (x). 
The Jacobian matrix of MD-GAD has the form : 



Ji 



1 

(j-2nnTj) -2(FTnI + nF^) 
L (J - 2nnTj - n^Jnl) 



(11) 



where, L is a n x n matrix involving the third order derivatives of V^. At any fixed point of the 
potential, since the derivative of V is zero, the term F^nl + nF^ is zero and the eigenvalues of Ji 
at these fixed points can be obtained from 



Ji 



1 

(j-2nnTj) 

L (J - 2nnTj - n^Jnl) 



(12) 
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This means, eigenvalues of Ji can be obtained from the eigenvalues of N = 
( J — n'^ Jnl — 2nn^ j) and square root of the eigenvalues of M = ( J — 2nn^ j) . 

Let Ai (x) < A2 (x) < A3 (x) < ... < (x) be the eigenvalues (including degeneracies) of 
H (x) and (wi, W2, w^) be the associated set of orthonormal eigenvectors at some fixed point 
of V . A negative (positive) eigenvalue of J (Hessian) corresponds to a stable manifold, while a 
positive eigenvalue of J corresponds to an unstable manifold. Assuming n equals wi (i.e. the 
minimum eigen mode of the local Hessian), we find that 

Mw, = (J - 2nn^j) w, = -A,w, + 2Ai5aw^ (13) 

Nw^ = (J - n'^Jnl - 2nn^j) -A^w^ + Aiw^ + 2\i5iiWi (14) 

Thus, the eigenvalues of Ji are — 2Ai} for i = 1 and —^/—Xi, (Ai — A^)} 

when _L wi (i.e. i > 1). At an index- 1 saddle point Ai < and A^ > for alH > 1. 
Consequently, the eigenvalues of Ji are all negative or complex numbers and this fixed point of 
V (x) becomes linearly stable attractor for the dynamical system 

At a locally stable fixed point of V (x) all the eigenvalues of H are positive. Thus, Ji has two 
positive eigenvalues, n negative eigenvalues and (2n — 2) complex eigenvalues. Correspondingly, 
there exists two unstable manifolds due to which this fixed point becomes unstable in MD-GAD. 
Detailed analysis of convergence of GAD to index- 1 saddle points for non-gradient systems has 
been reported elsewhere.'^ 

Now, for the dynamical system ^ the Jacobian near a fixed point of V (x) is given by 

r (J - 2nnTj) 
J2= ^ ^ (15) 

L (J - 2nnT J - Jnl) 

Following a similar procedure, the eigenvalues of the J2 are given by the eigenvalues of N and 
M. We find that all the eigenvalue of J 2 are negative at index- 1 saddle point of V (x). In contrast, 
there exists two positive and (2n-2) negative eigenvalues of J2 at a stable fixed point of V (x). 
Thus local minimum points of V (x) are linearly unstable while index- 1 saddle points become 
locally stable attractors. 
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VI. APPENDIX II:LOCAL CONVERGENCE ANALYSIS FOR ART 



Having analyzed the stability of different fixed points of V (x) within GAD, let us try to analyze 
the stability of fixed points on V (x) in ART. Originally, ART was proposed to explore the PES 
and only in subsequent modified versions it converges to a saddle pointP^^ Here we evaluate the 
propensity of ART to visit an index- 1 saddle point while exploring the PES. 

In ART, the system evolves by the following equation: 



X = F (x) - (F (x) , n) n, z/ > 1 



(16) 



where, n is the direction vector given by the difference of the current location x and initial location 
Xo in the configuration space, i.e. n = (x — Xq) / |x — Xo|. Close to a fixed point of V (x) since 



the force is zero, the Jacobian, J3, for (16) is J3 = (J — z/nn^j) where, J = — H = —V'^V (x). 
To illustrate the stability of the dynamical system in ([16]), let us consider for example a simple 2D 
potential energy landscape : 

V (x, y) = sin (ttx) sin (ny) . (17) 

At a saddle point, for instance at (0, 0), the Hessian has eigenvalues {1,-1}. Assuming n = 
[rii 712]^, with nl + nl = 1, the Jacobian of ART at a fixed point of V (x) is 



1 
1 



— jy 



nin2 



nin2 77^2 



At (0, 0), if J is given by 



-1 
1 



(1 — vn\) —unin2 
vnin2 (1 — yn\) 



(18) 



(19) 



Diagonalizing J3, we obtain the following eigenvalues 



(2-^) 



--J(l-8n?ni), 



^ + i;/(l-8n?ni) 



(20) 



Since, 1 = (nf + n|)^ > An\nl =^ n\nl < 1/4, consequently {1 — Snlnl) > —1. Thus, 
depending on the direction vector n the saddle can be stable or unstable for the dynamical equation 
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indlel). 
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FIG. 4: 4(a)| Shows the eigenvalues of the Hessian, arranged in ascending order, at the local minimum 



configuration for the vacancy ad-atom system. The upper inset shows that the Hessian is diagonal dominant 
and the lower inset shows the clustering of low lying eigenvalues of the Hessian. |4 (6) [ Shows the different 
DOF of the unit force vector (blue) and the superposition of 6 lowest eigenvectors of the Hessian at the local 
minimum configuration. The vulnerable DOF can be easily isolated. 
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FIG. 5: Simulation cell set-up for nanoindentation on (111) surface of an initially defect- free Cu thin-film 
containing 112,320 atoms. 5 (a) [ shows a cross section of the system after indentation performed using a 
spherical smooth indenter of radius 25 A on the (111) free surface. The indentation direction is perpendic- 
ular to the (111) surface. The atoms have been colored by their shear strain (dark blue corresponds to zero 
strain and red corresponds to maximum shear strain ^ 0.20). |5 (6) [ shows the contributions from different 
DOF to the unit force thus showing the atoms vulnerable to deformation. The inset shows a magnified 
portion of the plot. The vulnerable DOF can be easily identified. 
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FIG. 6: Some converged saddle configurations for tlie vacancy ad-atom system obtained from GAD sim- 
ulations. |6(a) I ad-atom diffusion (barrier 0.01 eV), |6 (b) | collective process involving ad-atom and vacancy 
migration (barrier 0.15 eV), |6(c)| vacancy ad-atom annihilation (barrier 0.30 eV), |6(d)l vacancy diffusion 
(barrier 0.59 eV), 6(e) [ sub-surface vacancy formation (barrier 0.77 eV), |6(/)] sub-surface vacancy forma- 
tion (barrier 0.82 eV). Some of the structures are shifted in plane for proper viewing. 
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FIG. 7: Some converged saddle configurations for the vacancy ad-atom system obtained from GAD sim- 
ulations. 7(a)| ad-atom formation by exchange mechanism (barrier 1.14 eV), |7(b)| divacancy and ad-atom 
formation (barrier 1.40 eV), |7(c)| divacancy and ad-atom formation (barrier 1.81 eV), \7{d)\ bulk- vacancy 
and ad-atom formation (barrier 2.08 eV). 
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FIG. 8: Summary of results depicting a wide spectrum of converged energy barriers for different rare events 
as a function of the distance of the final saddle configuration from the initial locally stable configuration for 
the surface vacancy and ad- atom system. The results correspond to simulations performed with qc = 0.50 
thus allowing the possibility of sampling a much larger portion of the configuration space for index- 1 saddle 
points. 
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FIG. 9: Summary of results depicting the distribution of converged saddle energies (with respect to the 
initial structure) for qc = 0.40 ( |9(a) ) and qc = 0.30 ( |9(6)| ) for the system with a surface vacancy and 
ad-atom. The results correspond to the evolution equations ([2]) 
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FIG. 1 1 : Some converged saddle configurations for the vacancy ad-atom configuration obtained from GAD 
simulations. |ll(a)| island sliding on (111) surface (barrier 0.39 eV), |l 1(b)] two atom sliding (barrier 0.49 
eV), |ll(c) single atom hopping (barrier 0.683), ll((j)| two atom sliding (barrier 0.69 eV), |l 1(e)] two atom 
sliding process taking place at two locations on the island (barrier 0.78 eV), |ll(/7 sliding of atoms (barrier 
0.84 eV). 
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(c) (d) 

FIG. 12: Some converged saddle configurations for the vacancy ad-atom configuration obtained from GAD 
simulations. 



12(a)| island rotation (barrier 0.92 eV), 12(5)| ad-atom formation on heptamer island (barrier 



1.39 eV), |12(c) two atom sliding process (barrier 1.49 eV), 12((j)| ad-atom and surface vacancy formation 



(barrier 1.86 eV). 
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FIG. 14: Summary of results depicting a wide spectrum of converged energy barriers for different rare events 
as a function of the distance of the final saddle configuration from the initial locally stable configuration for 
the system with a heptamer island on (111) surface. 
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